In our approach, first, the λ-hybrid potential is constructed by linearly combining the AA potential and CG potential as in Eq. (1). The λ satisfies and acts like a strength tuning parameter to reduce the ruggedness of the potential energy surface.

The free energy landscape of the hybrid system in the CG space, i.e., the PMF, is

To ensure the proper rigidity of bond and angle connection, the hybridization usually performs on part of AA potentials, such as the non-bonded term only, which includes the Van der Waals (VDW) potential and the Coulomb interaction, thus the total potential in the hybrid model is

A strongly hybrid potential (large λ) can reduce the atomic friction and speed up the simulation, while a lightly hybrid potential ( ) can provide more accurate atomic detail. A simple way to balance between the two points can communicate multiple λ-hybrid potentials involving λ=0. This communication can be implemented through the replica exchange technique.^[26] However, the hardware requirement of the replica exchange method is harsher than that of a conventional simulation, which will be n times of the conventional resources. Besides, it could encounter a problem in a low exchange rate. Here we apply another approach, the integrated temperature sampling (ITS),^[27] to linearly combine n hybrid potentials to a single effective potential,

In this paper, we illustrate the hybrid CG approach in proteins, where the average potential for each CG pair is fitted by the Lennard–Jones (LJ) 12–6 potential and the Coulomb interaction. To mimic the implicit solvent effect with ions, the Coulomb interaction is modified as screened Coulomb function with the relative dielectric constant , and the screen length l _D = 7 Å. In the current implementation of hyCG, instead of strictly parameterizing the average potential, we directly use the non-bonded parameter of the MARTINI force field^[29] as the CG potential. There are several features in MARTINI for supporting this tentative treatment. On the one hand, MARTINI aims at providing a general force field by assigning a phenomenological chemical function for CG beads. The attractive and the repulsive properties are modeled as the LJ potential, which is already fitted to including the VDW and the electric dipole interaction among various chemical environments. On the other hand, the mapping protocol is fine enough to describe the topology structure of the side chain, like the benzene ring in phenylalanine residue. Hence, it can avoid the overlap between different CG groups caused by atomistic diffusion while we largely reduce the atomistic force. We will adopt the latest version 2.2 of MARTINI for our model, which has been improved to correct certain hydrophobic interaction.^[30]

However, the closest distance of the LJ potential in MARTINI is a preassigned value, which in our test is found to be more repulsive for polar groups while we integrate it into our hybrid model. To avoid this repulsive force breaking the stability of a native globular protein, any pairs of the CG site within 4 Å in a native structure are excluded from the calculation of the LJ potential. The cut-off 4 Å is chosen to be slightly larger than the distance of two neighbor atoms so that the originally excluded list in MARTINI can be maintained and a very small fraction of non-bonded CG sites will be blind to each other. This treatment allows the energy to become further lower in the native structure.

The solvent effect in our simulation is treated as GB/SA implicit solvent^{[31, 32]} with surface tension setting to 0.005 kcal/mol/Ų and α-cutoff adopted 10 Å. The ion concentration in GB/SA is set as 0.2 M and all the simulations are carried out at 300 K and adopting charm22-cmap AA force field.^{[33, 34]} For the AA force field, the cut-off for both VDW and Coulomb is set as 12 Å and we apply a switched form for VDW. We adopt the same force field treatment for our average potential. The is chosen as 0.5.

The stability of protein’s native structure is critical for testing the accuracy of a newly developed force field. To understand whether the hyCG can stabilize native structures of proteins, a 35-residue Villin headpiece protein (PDB entry-code 1YRF) is used to test. We run a hyCG simulation and an AA simulation from the x-ray structure of 1YRF. Both are carried out for 100 ns.

The time evolution of the root of mean squared deviation (RMSD) of the protein backbone related to its native structure, shown in Fig. 1, indicates that the AA force field can stabilize the native structure up to 100 ns with the RMSD value below 3 Å. The stability of the hyCG force field can reach nearly 70 ns, and then the protein loses its tertiary structure but maintains a roughly correct secondary structure. The reduced stability of the hyCG indicates a lower potential barrier between the folded and the unfolded states. We select two pairs of distances, the end-to-end distance and the coil-to-coil distance (see Fig. 1(d)), to test the results of hyCG. The distributions of the distances in the hyCG and the AA (figs. 1(b) and 1(c)) are roughly overlapped with each other. The small deviation in the probability distribution of the hyCG in the CG space can be further refined by considering the difference between the hyCG effective potential in Eq. (7) and the AA potential based on reweighting.^[27]

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) Native stability of 1YRF tested by hyCG and AA. (a) Backbone RMSD, in which the red represents hyCG and green represents AA. The color is also adopted by the other figure. The black line indicates a 3 Å limit for the backbone RMSD threshold. Within this threshold, it is commonly considered as a folded state. (b) and (c) Distribution of two selected distances constructed from the first 50 ns trajectory. (d) Schematic representation of 1YRF native structure and the two selected distances.

To understand the formation of the secondary structure, we study the 1YRF protein again. We perform 250 ns simulation with the hyCG and the AA starting from a complete extended configuration. In order to understand how the non-bonded potential affects the formation efficiency of the secondary structure, we also carry out an energy enhanced sampling simulation on the non-bonded potential as a comparison, which is implemented by selective integrated temperature sampling (SITS).^[35]

Two collective variables are chosen to illustrate the formation rate of the helical structure and the correctness of the secondary structure, see Fig. 2. For both collective variables, the hyCG shares a similar growth in the secondary structure formation with the AA but is more efficient and stable. In comparison, although the SITS has a fast formation at the beginning, it will also visit many unfavorable regions for its ability in largely broadening the configuration. The overall effect will slow down the correct formation of the secondary structure in the SITS, which in turn indicates that the rate of the hyCG is most efficient. However, the percentage of the helix in the hyCG shows slightly over helical, yielding the most occupied state differing from the native structure during the short time simulations.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. (color online) Secondary structure formation of 1YRF using hyCG, SITS, and AA. All the secondary structure is determined by STRIDE in VMD.[46] (a) The percentage of helix structure. It is determined by the helical (STRIDE code H,G,I) number divided by the total residue number. The black line indicates the correct percentage of the helix number. (b) The correctness of the secondary structure. It is determined by the percentage of the correct secondary structure compared to the native structure.

To obtain enough samples for constructing the landscape in a short time simulation, we study a smaller protein instead of Villin. Here we study the folding of 20-residues Trp-cage protein (PDB entry-code 1L2Y) to understand the tertiary structure formation efficiency in the hyCG. Trp-cage is a widely studied protein due to its small size. To compare the efficiency of the hyCG and the AA force field, we first run N = 4 hyCG and AA simulations at 300 K, respectively, both starting at four different fully extended configurations generated at an equilibrium simulation under 700 K. Each trajectory has τ = 250 ns length. Among these, we find that there are n=2 hyCG trajectories reaching the native structure (with their backbone RMSD smaller than 2 Å) at the time and 90 ns, respectively, while the other 2 hyCG trajectories do not fold to the native structure within the τ-length simulations. Thus the folding time in the hyCG is estimated by^{[36, 37]}

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) Folding trajectory of 1L2Y with hyCG force field and AA force field. (a)–(d) Backbone RMSD of hyCG and AA. The initial configurations are all extended conformation taken from 700 K simulation. Here 2 Å threshold is enough to judge whether the protein is folded. (e) The best folding structure taken from the second trajectory of hyCG with the RMSD equals to 0.88 Å. The red one is the native structure and the blue one is the hyCG’s result.

To investigate whether the folding path is still preserved in the hyCG, we construct the metastable states and their corresponding transitional relationship. A folding path can be described as a transitional network^[39] from the unfolded state to the multiple metastable states and to the native state, which can be constructed by the current Markov state model^[40–42] or the trajectory mapping method.^{[43, 44]} Each metastable state and the main transitional path are visually shown in Fig. 5. We also exhibit two free energy surfaces (FES) to understand the roughness of the folding landscape. Each of the free energy surfaces is constructed by two collective variable sets: the fraction of native contact against the correctness of the secondary structure and the fraction of native contact against the backbone RMSD (reference is the native structure), respectively, see Fig. 4. The free energy is determined through the Markov state model.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. (color online) Folding landscape of 1L2Y along two given collective variables at 300 K. (a), (c) The fraction of native contact and the correctness of secondary structure. (d), (d) The fraction of native contact and the backbone RMSD (reference is the native structure). The energy is in units of kcal/mol. The free energies of both force fields are constructed by the Markov state model using multiple trajectories. The initial conformations are drawn from a short replica exchange simulation and it contains several nearly folded structures which can enhance the bridging between the folded and the intermediate states in normal simulation.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Schematically showing the folding path of 1L2Y. Each representative structure is identified by the Markov state model. The black arrows are marked to show the highly frequent transition events. Those states with small transitional rate and very short life are neglected in this diagram for the sake of simplicity.

Although some positions of FES local minimums in the hyCG fail to match strictly to that in the AA model, still, the main topology structure is maintained. Besides, the representative structure of important metastable states has a similar secondary and tertiary structure. Further, it verifies the similarity of their folding paths. However, certain metastable states, the V state in Fig. 5 with the backbone RMSD around 7 Å and the fraction of native contact between 0.75 and 0.80 in Fig. 4, are shown to be a branching folding path toward the native structure in the AA and now disappear in the hyCG. These metastable states contain a well-formed helix structure with a very loose coil-like structure, which will slow down the correct contact of the residue. Other simulations starting from this state with the hyCG indicate that the metastable structure in the AA is not stable anymore in the hyCG, see Fig. 6. In turn, it indicates a smoother folding landscape and a more efficient contact rate in the hyCG.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Stability of the branching metastable state (V state shown in Fig. 5) in the hyCG force field. The backbone RMSD is determined by comparing to the native structure. The black line and blue line are starting at two randomly selected conformations from the V metastable state. The red line represents the average backbone RMSD of the V state.